namespace Tuna {

template<class Tprec, int Dim>
inline bool Quick_Hay<Tprec, Dim>::calcCoefficients1D() 
{
    prec_t Gamma_dx = Gamma / dx;
    prec_t dx_dt = dx / dt;
    prec_t ce, cem, cep, cw, cwm, cwp, CE, CW;
    aE = 0.0; aW = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {
        CE = ce = u(i  ); 
	CW = cw = u(i-1); 
	cem = cep = 0.0;
	cwm = cwp = 0.0;

	// QUICK as presented in Hayase et al. J. of Comput. Phys., 98, 108-118, 1992.

// --- X
	if ( ce > 0 ) { 
	  CE = 0;
	  if (i == bi) {
	    cep = ce * (phi_0(i+1) - phi_0(i-1)) / 3.0;
	  } else {
	    cep = ce * 0.125 * (-phi_0(i-1) - 2*phi_0(i) + 3*phi_0(i+1));
	  }
	} else {
	  // The case i == ei is taken in to account in applyBoundaryConditions1D.
	  if (i == ei-1) {
	    cem = ce * (phi_0(i+2) - phi_0(i)) / 3.0;
	  } else if (i < ei-1) {
	    cem = ce * 0.125 * (-phi_0(i+2) - 2*phi_0(i+1) + 3*phi_0(i));
	  }
	}

	if ( cw > 0 ) {
	  // The case i == bi is taken in to account in applyBoundaryConditions1D.
	  if (i == bi+1) {
	    cwp = cw * (phi_0(i) - phi_0(i-2)) / 3.0;
	  } else if (i > bi+1) {
	    cwp = cw * 0.125 * (-phi_0(i-2) - 2*phi_0(i-1) + 3*phi_0(i));
	  }
	} else {
	  CW = 0;
	  if (i == ei) {
	    cwm = cw * (phi_0(i-1) - phi_0(i+1)) / 3.0;
	  } else {
	    cwm = cw * 0.125 * (-phi_0(i+1) - 2*phi_0(i) + 3*phi_0(i-1));
	  }
	}

	aE (i) = Gamma_dx - CE;
	aW (i) = Gamma_dx + CW;
	aP (i) = aE (i) + aW (i) + dx_dt + (ce - cw) ;
	    
	sp (i) = phi_0(i) * dx_dt - (cep + cem - cwp - cwm) ;
    }

    /*** *
    // Dirichlet Boundary Condition on LEFT WALL
    prec_t phi_A = 1;
    aP(bi) += aW(bi) + Gamma_dx;
    aE(bi) += Gamma_dx / 3.0;  
    sp(bi) += (2 * aW(bi) + 2.0 * Gamma_dx / 3.0) * phi_A;
    aW(bi) = 0.0;

    /*** /
    // Dirichlet Boundary Condition on LEFT WALL
    prec_t phi_B = 0;
    aP(ei) += aE(ei) + Gamma_dx;
    aW(ei) += Gamma_dx / 3.0;  
    sp(ei) += (2 * aE(ei) + 2.0 * Gamma_dx / 3.0) * phi_B;
    aE(ei) = 0.0;

    /* ****/

    applyBoundaryConditions1D();
    return 0;
}
          
template<class Tprec, int Dim>
inline bool Quick_Hay<Tprec, Dim>::calcCoefficients2D() 
{
    prec_t Gdy_dx = Gamma * dy / dx;
    prec_t Gdx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    prec_t ce, cem, cep, cw, cwm, cwp, CE, CW;
    prec_t cn, cnm, cnp, cs, csm, csp, CN, CS;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;



    for (int i =  bi; i <= ei; ++i)
      for (int j = bj; j <= ej; ++j)
	{
	  CE = ce = u(i, j) * dy;
	  CW = cw = u(i-1, j) * dy;
	  CN = cn = v(i, j) * dx;
	  CS = cs = v(i, j-1) * dx;
	  cem = cep = 0.0;
	  cwm = cwp = 0.0;
	  cnm = cnp = 0.0;
	  csm = csp = 0.0;

	  // QUICK as presented in Hayase et al. J. of Comput. Phys., 98, 108-118, 1992.
// --- X
	  if ( ce > 0 ) { 
	    CE = 0;
	    if (i == bi) {
	      cep = ce * (phi_0(i+1,j) - phi_0(i-1,j)) / 3.0;
	    } else {
	      cep = ce * 0.125 * (-phi_0(i-1,j) - 2*phi_0(i,j) + 3*phi_0(i+1,j));
	    }
	  } else {
	  // The case i == ei is taken in to account in applyBoundaryConditions2D.
	    if (i == ei-1) {
	      cem = ce * (phi_0(i+2,j) - phi_0(i,j)) / 3.0;
	    } else if (i < ei-1){
	      cem = ce * 0.125 * (-phi_0(i+2,j) - 2*phi_0(i+1,j) + 3*phi_0(i,j));
	    }
	  }
	  
	  if ( cw > 0 ) {
 	  // The case i == bi is taken in to account in applyBoundaryConditions2D.
	    if (i == bi+1) {
	      cwp = cw * (phi_0(i,j) - phi_0(i-2,j)) / 3.0;
	    } else if (i > bi+1) {
	      cwp = cw * 0.125 * (-phi_0(i-2,j) - 2*phi_0(i-1,j) + 3*phi_0(i,j));
	    }
	  } else {
	    CW = 0;
	    if (i == ei) {
	      cwm = cw * (phi_0(i-1,j) - phi_0(i+1,j)) / 3.0;
	    } else {
	      cwm = cw * 0.125 * (-phi_0(i+1,j) - 2*phi_0(i,j) + 3*phi_0(i-1,j));
	    }
	  }
// --- Y
	  if ( cn > 0 ) { 
	    CN = 0;
	    if (j == bj) {
	      cnp = cn * (phi_0(i,j+1) - phi_0(i,j-1)) / 3.0;
	    } else {
	      cnp = cn * 0.125 * (-phi_0(i,j-1) - 2*phi_0(i,j) + 3*phi_0(i,j+1));
	    }
	  } else {
	    if (j == ej-1) {
	      cnm = cn * (phi_0(i,j+2) - phi_0(i,j)) / 3.0;
	    } else if (j < ej-1) {
	      cnm = cn * 0.125 * (-phi_0(i,j+2) - 2*phi_0(i,j+1) + 3*phi_0(i,j));
	    }
	  }
	  
	  if ( cs > 0 ) { 
	    if (j == bj+1) {
	      csp = cs * (phi_0(i,j) - phi_0(i,j-2)) / 3.0;
	    } else if (j > bj+1) {
	      csp = cs * 0.125 * (-phi_0(i,j-2) - 2*phi_0(i,j-1) + 3*phi_0(i,j));
	    }
	  } else {
	    CS = 0;
	    if (j == ej) {
	      csm = cs * (phi_0(i,j-1) - phi_0(i,j+1)) / 3.0;
	    } else {
	      csm = cs * 0.125 * (-phi_0(i,j+1) - 2*phi_0(i,j) + 3*phi_0(i,j-1));
	    }
	  }

	  aE (i,j) = Gdy_dx - CE;
	  aW (i,j) = Gdy_dx + CW;
	  aN (i,j) = Gdx_dy - CN;
	  aS (i,j) = Gdx_dy + CS;
	  aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt
	    + (ce - cw) + (cn - cs) ;	    
	  sp (i,j) = phi_0(i,j) * dxy_dt
	    - (cep + cem - cwp - cwm + cnp + cnm - csp - csm);	    
	}    
    applyBoundaryConditions2D();
    return 0;
}


//
//---------------------  3D  ---------------------
//
template<class Tprec, int Dim>
inline bool Quick_Hay<Tprec, Dim>::calcCoefficients3D()
{
    prec_t dyz = dy * dz, dyz_dx = Gamma * dyz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dxz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dxy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    prec_t ce, cem, cep, cw, cwm, cwp, CE, CW;
    prec_t cn, cnm, cnp, cs, csm, csp, CN, CS;
    prec_t cf, cfm, cfp, cb, cbm, cbp, CF, CB;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int k = bk; k <= ek; ++k)
      for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	  {
	    CE = ce = u(i  , j, k) * dyz;
	    CW = cw = u(i-1, j, k) * dyz;
	    CN = cn = v(i, j  , k) * dxz;
	    CS = cs = v(i, j-1, k) * dxz;
	    CF = cf = w(i, j, k  ) * dxy;
	    CB = cb = w(i, j, k-1) * dxy;
	    cem = cep = 0.0;
	    cwm = cwp = 0.0;
	    cnm = cnp = 0.0;
	    csm = csp = 0.0;
	    cfm = cfp = 0.0;
	    cbm = cbp = 0.0;

	    // QUICK as presented in Hayase et al.
// --- X
	    if ( ce > 0 ) { 
	      CE = 0;
	      if (i == bi) {
		cep = ce * (phi_0(i+1,j,k) - phi_0(i-1,j,k)) / 3.0;
	      } else {
		cep = ce * 0.125 * (-phi_0(i-1,j,k) - 2*phi_0(i,j,k) + 3*phi_0(i+1,j,k));
	      }
	    } else {
	      // The case i == ei is taken in to account in applyBoundaryConditions3D.
	      if (i == ei-1) {
		cem = ce * (phi_0(i+2,j,k) - phi_0(i,j,k)) / 3.0;
	      } else if (i < ei-1) {
		cem = ce * 0.125 * (-phi_0(i+2,j,k) - 2*phi_0(i+1,j,k) + 3*phi_0(i,j,k));
	      }
	    }
	    
	    if ( cw > 0 ) { 
	      // The case i == bi is taken in to account in applyBoundaryConditions3D.
	      if (i == bi+1) {
		cwp = cw * (phi_0(i,j,k) - phi_0(i-2,j,k)) / 3.0;
	      } else if (i > bi+1) {
		cwp = cw * 0.125 * (-phi_0(i-2,j,k) - 2*phi_0(i-1,j,k) + 3*phi_0(i,j,k));
	      }
	    } else {
	      CW = 0;
	      if (i == ei) {
		cwm = cw * (phi_0(i-1,j,k) - phi_0(i+1,j,k)) / 3.0;
	      } else {
		cwm = cw * 0.125 * (-phi_0(i+1,j,k) - 2*phi_0(i,j,k) + 3*phi_0(i-1,j,k));
	      }
	    }

// --- Y
	    if ( cn > 0 ) { 
	      CN = 0;
	      if (j == bj) {
		cnp = cn * (phi_0(i,j+1,k) - phi_0(i,j-1,k)) / 3.0;
	      } else {
		cnp = cn * 0.125 * (-phi_0(i,j-1,k) - 2*phi_0(i,j,k) + 3*phi_0(i,j+1,k));
	      }
	    } else {
	      if (j == ej-1) {
		cnm = cn * (phi_0(i,j+2,k) - phi_0(i,j,k)) / 3.0;
	      } else if (i < ei-1) {
		cnm = cn * 0.125 * (-phi_0(i,j+2,k) - 2*phi_0(i,j+1,k) + 3*phi_0(i,j,k));
	      }
	    }
	    
	    if ( cs > 0 ) { 
	      if (j == bj+1) {
		csp = cs * (phi_0(i,j,k) - phi_0(i,j-2,k)) / 3.0;
	      } else if (j > bj+1) {
		csp = cs * 0.125 * (-phi_0(i,j-2,k) - 2*phi_0(i,j-1,k) + 3*phi_0(i,j,k));
	      }
	    } else {
	      CS = 0;
	      if (j == ej) {
		csm = cs * (phi_0(i,j-1,k) - phi_0(i,j+1,k)) / 3.0;
	      } else {
		csm = cs * 0.125 * (-phi_0(i,j+1,k) - 2*phi_0(i,j,k) + 3*phi_0(i,j-1,k));
	      }
	    }

// --- Z
	    if ( cf > 0 ) { 
	      CF = 0;
	      if (k == bk) {
		cfp = cf * (phi_0(i,j,k+1) - phi_0(i,j,k-1)) / 3.0;
	      } else {
		cfp = cf * 0.125 * (-phi_0(i,j,k-1) - 2*phi_0(i,j,k) + 3*phi_0(i,j,k+1));
	      }
	    } else {
	      if (k == ek-1) {
		cfm = cf * (phi_0(i,j,k+2) - phi_0(i,j,k)) / 3.0;
	      } else if (k < ek-1) {
		cfm = cf * 0.125 * (-phi_0(i,j,k+2) - 2*phi_0(i,j,k+1) + 3*phi_0(i,j,k));
	      }
	    }
	    
	    if ( cb > 0 ) { 
	      if (k == bk+1) {
		cbp = cb * (phi_0(i,j,k) - phi_0(i,j,k-2)) / 3.0;
	      } else if (i > bk+1) {
		cbp = cb * 0.125 * (-phi_0(i,j,k-2) - 2*phi_0(i,j,k-1) + 3*phi_0(i,j,k));
	      }
	    } else {
	      CB = 0;
	      if (k == ek) {
		cbm = cb * (phi_0(i,j,k-1) - phi_0(i,j,k+1)) / 3.0;
	      } else {
		cbm = cb * 0.125 * (-phi_0(i,j,k+1) - 2*phi_0(i,j,k) + 3*phi_0(i,k,k-1));
	      }
	    }

	    aE (i,j,k) = dyz_dx - CE;
	    aW (i,j,k) = dyz_dx + CW;
	    aN (i,j,k) = dxz_dy - CN;
	    aS (i,j,k) = dxz_dy + CS;
	    aF (i,j,k) = dxy_dz - CF;
	    aB (i,j,k) = dxy_dz + CB;
	    aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + aN (i,j,k) + aS (i,j,k)
	      + aF (i,j,k) + aB (i,j,k) + dxyz_dt
	      + (ce - cw) + (cn - cs) + (cf - cb);

	    sp (i,j,k) = phi_0(i,j,k) * dxyz_dt
	      - (cep + cem - cwp - cwm + 
		 cnp + cnm - csp - csm +
		 cfp + cfm - cbp - cbm);	  
	}    
    applyBoundaryConditions3D();
    return 0;
}

} // Tuna namespace
















